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ABSTRACT 



A method for rapidly forming a stochastic model of Gaussian or related type, 

10 representative of a porous heterogeneous medium such as an underground reservoir, 
constrained by data characteristic of the displacement of fluids. The method comprises 
construction of a chain of realizations representative of a stochastic model (Y) by 
gradually combining an initial realization of (Y) and one or more other realization(s) of 
(Y) referred to as a composite realization, and minimizing an objective function (J) 

15 measuring the difference between a set of non-linear data deduced from the 
combination by means of a simulator simulating the flow in the medium and the data 
measured in the mediimi, by adjustment of the coefficients of the combination. The 
composite realization results from the projection of the direction of descent of the 
objective function, calculated by the flow simulator for the initial realization, in the 

20 vector subspace generated by P realizations of (Y), randomly drawn and independent of 
one another, and of the initial realization. During optimization, the chain is explored so 
as to identify a realization that allows minimizing the objective function (J). In order to 
sufficiently reduce the objective function, sequentially constructed chains are explored 
by taking as the initial realization the optimum realization determined for the previous 

25 chain. The method may be used for development of oil reservoirs . 
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BACKGROUND OF THE INVENTION 
Field of the Invention 

The present invention relates to a method for more rapidly forming a stochastic 
numerical model of Gaussian or a related type, representative of a porous heterogeneous 
5 medium (such as a hydrocarbon reservoir for example) calibrated in relation to data 
referred to as dynamic data, characteristic of the displacement of fluids in the medium 
such as, for example, production data (pressures obtained from well tests, flow rates, 
etc.). 



1 0 Description of the Prior Art 

Optimization in a stochastic context consists in determining realizations of a 
stochastic model which meet a set of data observed in the field. In reservoir 
engineering, the realizations to be identified correspond to representations, in the 
reservoir field, of the distribution of carrying properties such as the permeability or 

15 porosity. These realizations form numerical reservoir models. The available data are, for 
example, isolated permeability or porosity measurements, a spatial variability model 
determined according to isolated measurements or data directly related to the fluid 
flows in an underground reservoir, that is pressures, breakthrough times, flow rates, etc. 
The data are often not linearly related to the physical properties to be modelled. A 

20 randomly drawn realization is generally not in accordance with the whole of the data 
collected. Coherence in relation to the data is integrated in the model by means of an 
inverse procedure See 
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- Tarantola, A., "Inverse Problem Theory - Methods for Data Fitting and Model 
Parameter Estimation", Elsevier Science Publishers, 1987. 

The simplest technique is therefore the trial-and-error method. This approach 
randomly takes realizations until a realization meeting the data is obtained. It affords the 
5 advantage of conservation of the spatial variability of the model but requires prohibitive 
calculation time. It is therefore very rarely used in practice. 

An option that is often preferred is based on gradients calculation. The gradients 
methods allow modifying an initial realization in a direction of search which is 
estimated from the gradients. The modifications are applied iteratively until data 
10 calibration is considered to be acceptable. The gradients methods are attractive because 
of their efficiency. However, they are no longer suitable when the number of 
parameters, that is the number of permeability and porosity values forming the 
numerical model, is large. Besides, they do not allow modifying the realizations while 
respecting the spatial structure of the stochastic model. 

15 More recently, a geostatistical parameterization technique has been introduced to 

constrain, by gradual deformation, the stochastic realizations to data on which they 
depend non-linearly. This technique is described in French patents 2,780,798 and 
2,795,841 of the assignee, and in the following publications, notably : 

- Hu, L.Y., 2000, Gradual Deformation and Interative Calibration of Gaussian-Related 
20 Stochastic Models : Math. Geology, Vol.32, No.l, 

- Le Ravalec, M. Et al., 2000, The FFT Moving Average (FFT-MA)Generator: An 
Efficient Numercial Method for Generating and Conditioning Gaussian Simulations: 
Math. Geology, Vol.32, No.6, 
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- Hu, L.Y., Blanc, G. And Noetinger, B. (2001) : Gradual Deformation and Interative 
Calibration of Sequential Stochastic Simulations. Math. Geology, Vol.33, No.4. 

This method has been successfully applied in various cases, notably from data 
obtained from oil development fields, as described in the following documents : 

5 - Roggero, F. et al., 1998, Gradual Deformation of Continuous Geostatistical Models for 
History Matching, paper SPE 49004 : Proc. SPE Annual Technical Conference and 
Exhibition, New Orleans, 

- Hu, L.Y. et al., 1998, Constraining a Reservoir Facies Model to Dynamic Data Using a 
Gradual Deformation Method, paper B-01 : Proc. 6th European Conference on 

10 Mathematics of Oil Recovery (ECMOR VI), 8-1 1 September 1998, Peebles, Scotland. 

As described in detail hereafter, the gradual deformation method allows gradual 
modification a realization of a stochastic model of Gaussian or a Gaussian-related type 
while respecting the spatial structure of the model. 

Stochastic optimization 

15 Let /^^^ = fi'^^-^ /m^'O be the field data and / = (/„ /2,..., /m) the 

corresponding responses simulated for a realization z of a given stochastic model Z. In 
general, the responses / = (/i, /a,. •» /m) are obtained by solving numerically the direct 
problem. Thus, if z represents a permeability field, data / can be pressure 
measurements. In this case, they are simulated from a flow simulator. The goal of a 

20 stochastic optimization is to produce realizations of Z which reduce the differences 
betweien the observed data and the numerically simulated corresponding responses. 
These differences are measured by the following objective function : 
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^ m=l 

Coefficients co^ are weights assigned to data f^^^, which are functions of 
realization z. In this sense, minimization of the objective function is a problem with 
several variables. 

5 Let N be the number of grid cells or of components forming realization z. N is often 

very large (10"^ ~ 10^). It is therefore very difficult to perform an optimization directly in 
relation to the components of z. Furthermore, realization z, even modified, must remain 
a realization of Z. The gradual deformation method allows these difficulties to be 
overcome. 

10 Random search from the gradual deformation method 

A random field Z is now considered that can be transformed into a Gaussian 
random field Y. The gradual deformation technique allows construction of a continuous 
chain of realizations by combining an initial realization yo of Y with another realization 
ui, referred to as complementary, of Y, ui being independent of yo (Figure la). The 
15 combination coefficients are for example cos(t) and sin(t), and the combined realization 
meets the relation : 

yW = yo cos t + ui sin t 
where t is the deformation parameter. 

Another realization chain construction technique combines the initial realization no 
20 longer with one, but with P complementary realizations Up (p=l,P) of Y (Figure lb). 
The coefficients of the combination are such that the sum of their squares is 1 . 
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Once the chain is formed, it can be explored by varying the deformation parameter t 
and while trying to identify, from among all the realizations of this chain, the realization 
which minimizes the objective function. This minimization is performed in relation to t. 
Parameterization according to the gradual deformation method allows reduction of the 
5 number of dimensions of the problem from N to 1, where N is the number of values that 
constitute the field to be constrained. Furthermore, the sum of the combination 
coefficients squared being 1, the optimized realization y(topt) still is a realization of Y : 
The optimized relation follows the same spatial variability model as all the realizations 
ofY. 

10 However, if the exploration of the realization space is restricted to a single chain, 

the possibilities of sufficiently reducing the objective function are greatly limited. The 
above procedure therefore has to be repeated, but with new realization chains. These 
realization chains are constructed successively by combining an initial realization which 
is here the optimum realization determined at the previous iteration, with a 

15 complementary realization of Y, randomly drawn each time. Thus, at iteration /, the 
continuous realization chain is written as follows : 

yi (t) = yi-i cos t + ui sin t . 

yM is the optimum realization defined at iteration /-I and the ui are independent 
realizations of Y. The latter are also independent of yo. This formulation implies that the 
20 realization chain corresponds to a hyperellipse of dimension N. 

Minimizing the objective function in relation to t allows improving or at least to 
preserve calibration of the data each time a new realization chain is explored. This 
iterative minimum search procedure is continued as long as data calibration is not 
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satisfactory. The random side of the method lies in the fact that, upon each iteration, the 
complementary realization is drawn at random. In fact, the direction of search that is 
followed from the optimized realization at the previous stage is random. The direction 
of search, for a given chain and from the optimum realization defined above, is : 



This direction of search only depends on ul. Furthermore, ul being independent of 
the complementary realizations already generated ui, U2,..., ui_i and also of yo, the 
direction of search at the start of each new chain is orthogonal to the tangent defined for 
the previous chain at the same point (Figure 2). Although it may seem appropriate to 
10 select a direction of search orthogonal to this tangent, there is an infinite number of 
possible orthogonal directions. 

Experience shows that, after several iterations, the new directions of search no 
longer contribute significantly to the decrease of the objective fiinction (Figure 6). 

It has also been considered to combine the initial realization not only with one, but 
15 with several complementary realizations. In this case, the number of deformation 
parameters increases and it is equal to ttie number of complementary realizations 
involved in a gradual combination. Although the optimization process is then more 
flexible, several parameters have to be managed, which is not easy. Besides, such a 
process is not necessarily more efficient because it can depend on the execution of a 
20 larger number of direct flow simulations. 




= -J^z-i sin 0 + cosO 
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SUMMARY OF THE INVENTION 

The mefhod according to the invention allows more rapid formation of a numerical 
model representative of the distribution of a physical quantity in a porous heterogeneous 
medium such as an underground zone (oil reservoirs, aquifers, etc.), constrained by data 
5 collected in the medium (dynamic data characteristic of the displacement of fluids in the 
medium), collected by measurements (in production, injection or observation wells for 
example) or previous observations. 

The method according to the invention has applications for underground zone 
modelling that generates representations showing how a physical quantity is distributed 
10 in a subsoil zone (permeability notably), best compatible with observed or measured 
data, in order for example to favor the development thereof 

It comprises an iterative process of gradual deformation wherein an initial 
realization of at least part of the selected model of the medium is linearly combined 
with at least a second realization independent of the initial realization, the coefficients 
15 of this linear combination being such that the sum of their squares is 1, and an objective 
function measuring the difference between a set of non-linear data deduced from the 
combination by means of a medium simulator and the geologic and dynamic data is 
minimized by adjusting the coefficients of the combination, the iterative process being 
repeated until an optimum realization of the stochastic model is obtained. 

20 With the method the rate of gradual deformation to the optimum model 

representative of the medium is accelerated by selecting as the second realization to be 
combined with the initial realization at least one composite realization obtained by 
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selecting beforehand a direction of descent defined as a function of the gradients of the 
objective function in relation to all the components of the initial realization. 

The composite realization is obtained for example by linear combination of a set of 
independent realizations of the model, the coefficients of the combination being 
5 calculated so that the direction of descent from the initial realization y is as close as 
possible to the realization defined by the gradients of the objective function in relation 
to all the components of the initial realization. 

Optimization is for example carried out from a deformation parameter which 
controls the combination between the initial realization and the composite realization. 

10 In cases where tfie combination affects orily part of the initial realization, the 

iterative process of gradual deformation is applied to Gaussian white noise used to 
generate a Gaussian realization and the derivatives of the objective function with 
respect to the components of the Gaussian white noise are determined. 

According to an implementation mode, the initial realization is combined with a 
15 certain number M of composite realizations, all obtained by composition from Pm 
independent realizations of Y, the optimization involving M parameters. 

In other words, the method comprises a new gradual combination scheme which 
accounts for the information provided by the gradients at the initial point of any chain of 
realizations. Construction of a chain always starts from an initial realization and from a 
20 set of complementary realizations, all independent and coming fi-om the same stochastic 
model. The initial realization is however not directly combined with the complementary 
realizations. The complementary realizations make it possible to explore the realizations 
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space in different directions. These directions are not equivalent : some better approach 
to the optimum. At this stage, a realization referred to as composite realization is 
elaborated by combining the complementary realizations only. A chain of realizations is 
then created from the initial realization and from this composite realization. This chain, 
like the chain proposed in the basic case of gradual deformation, can be explored by 
means of a single deformation parameter. 

The composite realization is constructed to provide a direction of search along the 
chain as close as possible to the direction of descent given by the gradients. As 
mentioned above, all complementary realizations are not equivalent : the composite 
realization takes the best of each complementary realization. 

The method allows reaching the formation of the numerical model representative of 
the medium more rapidly. 

BRIEF DESCRIPTION OF THE FIGURES 

Other features and advantages of the method according to the invention will be 
clear from reading the description hereafter of an application given by way of non 
limitative example, with reference to the accompanying drawings wherein : 

- Figures la, lb show gradual deformation schemes (referred to as GDM) that are 
already known, 

- Figure Ic shows the gradual deformation scheine (GBC) corresponding to the method 
according to the invention, 

- Figure 2 shows realization chains in a Euclidean space with N dimensions, where the 
tangent at the level of the optimized realization for a realization chain /-I (RC/-i) is 



11 

orthogonal to the direction of search for the initial realization of the next realization 
chain / (RCi), 

- Figure 3 shows the projection of the direction of descent v in the subspace U defined 
by the orthonormal base formed by P independent realizations (ui,..., Up), 

5 - Figures 4A to 4C show permeability distribution examples respectively for the 
reference realization, the initial realization and the realization constrained to the 
pressure data, 

- Figures 5A to 5E respectively show the variations as a function of time of the 
bottomhole pressures respectively simulated for the five wells (BHP-OBSl, BHP- 

10 OBS2, BHP-PROl, BHP-OBS3, BHP-OBS4) of Figures 4, respectively for the 
reference (data), initial and constrained (match) permeability distributions, and 

- Figure 6 shows the evolution (OF) of the objective function as a function of the 
number k of flow simulations performed, GDMl corresponding to an optimization 
carried out by combining an initial realization and a single complementaiy realization, 

15 GDMGBC3, to an optimization carried out by combining the initial realization and a 
composite realization constructed from three complementary realizations, GDMGBCIO, 
to the optimization carried out by combining the initial realization and a composite 
realization constructed from ten complementary realizations, and GDMGBC30, to an 
optimization carried out by combining the initial realization and a composite realization 

20 constructed from thirty complementary realizations. 

DETAILED DESCRIPTION OF THE INVENTION 

The method according to the invention allows orienting each iteration of the 
construction of the realization chain in order to reach a desirable direction of search. 
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The technique selected takes advantage of the information provided by the direction of 
descent defined by the gradients of the objective function J. The technique can be 
implemented by means of a numerical sihiulator of a type known to the man skilled in 
the art, such as the ATHOS or ECLIPSE simulators. 

Search oriented by the direction of descent 

At this stage consideration is given to the direction of descent (evolution towards a 
minimum value, whether local or not) defined by the gradients of the objective fimction 
J in relation to all the components of realization z. These gradients are deduced firom 
sensitivity coefficients 9/m / 9zn : 

The problem of sensitivity coefficients calculation has been widely dealt with in the 
scientific literature. The following document can for example be referred to : 

- Sun, N.Z., Inverse Problems in Groundwater Modelling, Kluwer Acad. Publ., 
Dordrecht, The Netherlands, 1994. 

In particular, the adjoint state technique allows calculation of all of coefficients by 
solving the direct problem and the adjoint problem, within a time interval equivalent to 
twice the time required for solution of the direct problem. 

Since the goal is to integrate the information provided by the gradients in the 
gradual deformation method, going back to the Gaussian realization is necessary which 
results from the transformation of z. Two cases can be distinguished according to 
whether isolated conditioning data, that is measurements of z at certain points, are 
available or not. 
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When no data are available, if y is the Gaussian realization obtained by 
transforming z, then : 

dJ ^ dJ dz„ 

In the opposite case, the assumption is that z is known at certain points a. Thus, the 
5 Gaussian realization s, deduced from the transformation of z, is a conditional realization 
of Y. To obtain s, a non-conditional realization y of Y is generated and conditioned to 
the values known at a by kriging. The conditional realization is deduced from : 

s = s* -H (y - y*). 

s* and y* are respectively calculated by kriging from the real data and the data of y 
10 generated at the level of points a. It can then be shown that : 

f dJ dz„ _ 
dJ -yVn^a 

[0, V« = a 

For continuous physical properties, dzn/dyn or dzn/dsn are calculated from the 
anamorphosis frmction allowing to transform realization z into a Gaussian realization. 
When the physical properties considered are category or discrete properties, these 
15 derivatives do not exist. The gradients techniques then cannot be applied. 

These various relations are of direct interest if a realization just has to be deformed 
globally. On the other hand, if it is desired to deform the realization locally, the gradual 
deformation method has to be applied to the Gaussian white noise x used to generate y. 
In this case, an additohal stage is required : calculation of the derivatives of the 
20 objective function with respect to the components of the Gaussian white noise. 
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To illustrate the calculation of these derivatives, concentration is made on the 
particular case where the non-conditional Gaussian realization y is obtained from the 
FFT-MA generator described in the article published by Le Ravalec et al. 2000 
mentioned above. 

The basic principle of this FFT-MA (FFT-Moving Average) generator is to 
transform a Gaussian white noise x into a Gaussian realization y correlated from a 
convolution product : 

y = g * X. 

Function g results from the expansion of covariance function C such that C = g * g\ 
where g' is the transpose of g. The derivatives of the objective fimction with respect to 
the components of the Gaussian white noise are : 

dJ _\ry dJ dy^ 

The discrete expression for the convolution product leads to * , which 

implies 9yi/9xn = gi-n. If this formulation is introduced in the derivatives of the 
objective function, it can be shown that : 



^ = V — 



This formula expresses the fact that the derivative of the objective function with 
respect to the n* component of the Gaussian white noise is given by the n* component 
of the field obtained by convoluting all the derivatives of the objective function with 
respect to the components of the Gaussian realization with the kernel of the covariance 
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function. From the framework established for the FFT-MA generator, these derivatives 
are determined as follows. 

1 -Calculation of the Fourier transform of dJ I 9y, these derivatives are obtained by 
means of the direct numerical simulator ; 

5 2-Multiplication of this Fourier transform with that of g that is provided by FFT- 

MA ; 

3-Calculation of the inverse Fourier transform of the previous product. 

The time required for calculation of these derivatives is negligible : it represents an 
additional time of 2/3 in relation to the simulation of a Gaussian realization by FFT- 
10 MA. 

Whatever the realization generator and the direct numerical simulator, it is assumed 
hereafter that a direction of descent from the derivatives of the objective frinction can be 
defined. If optimization of the objective function is performed in relation to this 
direction of descent only, the coherence of the realization in relation to the spatial 
15 variability model is generally destroyed. In the section hereunder, the information 
provided by the derivatives of the objective function is integrated in the gradual 
combination scheme. 

Taking account of the derivatives of the objective function in the gradual 
deformation process 

20 The realization chain yl(t) constructed from yO and from another realization u of Y 

(Figure la) is considered. Now, instead of using a complementary realization u as it is, 
P complementary realizations Up (p=l,2,...,P) of Y are randomly drawn and u is 
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substituted for a combination of the Up=i,p (Figure Ic). This combination is not any 
combination : it follows the following construction mode : 

p p 

The resulting realization u is a realization of Y and it is independent of yO. 
According to Eq.l, u is also the direction of search calculated for chain yi(t) at the 
starting point yo. The method constructs u so that the direction u is as close as possible 
to the direction of descent given by the gradients at yo. 

The space V of the vectors is first defined with N dimensions provided with the 
scalar product : 

yi n and yj.„ are respectively the n* components of vectors yj and yjU is defined as a 
subspace of V defined by the orthonormal base (ui,U2,...,Up). The direction of search in U 
which is the closest to the direction of descent v is given by the projection of v in U 
(Figure 3) : 

p 

By normalizing this vector, the desired direction u is obtained. The combination 
coefficients \ of Eq.(2) are thus : 

The realization u thus defined is referred to as composite realization. 

So far the construction of realization chains by combination of the initial realization 
with a composite realization has been considered. This technique can however be 
generalized to the construction of chains involving a certain number M of composite 
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ZM p 
. , J w.w^w 

realizations being independent. Optimization then involves M parameters. This 
technique increases the degree of freedom in the optimization process, but requires 
management of M optimization parameters. 

5 Numerical example 

A synthetic reservoir model is constructed on which the method according to the 
invention is tested. 

The synthetic reference reservoir is shown in Figure 4A. It is a monolayer reservoir 
comprising 51x51 grid cells which are 10 m thick and 40 m in side. The permeability 

10 distribution is lognormal with an average of 200 mD and a standard deviation of 100 
mD. The logarithm of the permeability field is characterized by an isotropic spherical 
variogram and a correlation length of 480 m. The other petrophysical properties are 
constant : the porosity is 25 %, the total compressibility 5.10"^ bar"* and the fluid 
viscosity 1 cP. A production well BHP-PROl, with a radius of 7.85 cm, fi:ee of any skin 

15 effect, is at the center of the reservoir and is surrounded by four observation wells 
(BHP-OBSl, BHP-OBS2, BHP-OBS3, BHP-OBS4) (Figures 4). A numerical flow 
simulation allows obtaining, for this reservoir, a set of reference data comprising the 
bottomhole pressures of the five wells (Figure 5). 

The object of the inverse problem is to determine a reservoir model coherent with 
20 the pressure data, the permeabilities distribution being assumed to be unknown. Four 
optimization processes are therefore launched, starting firom the same initial realization 
(Figure 4B). For each process, a single optimization parameter is considered, that is the 
deformation parameter. The first process (GDMl) takes up the conventional gradual 
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deformation scheme with construction of a realization chain using the initial realization 
and a single complementary realization. The other three processes (GDMGBC3, 
GDMGBCIO and GDMGBC30) illustrate the application of the method according to 
the invention : the chains are in this case elaborated with the initial reaUzation and a 

5 composite realization obtained from the combination of 3, 10 and 30 complementary 
realizations. The composite realization is constructed as explained above (Eq.2) by 
integrating the information provided by the gradients of the objective function with 
respect to the Gaussian white noise. For each process, the evolution of the objective 
function (OF) as a function of the number k of flow simulations performed is shown in 

10 Figure 6. 

It can be observed that using the gradients and increasing the number of 
complementary realizations allows the objective function to decrease more rapidly for 
the same number of flow simulations. 



